! ----------------------------
!
! Main program.
!
! ----------------------------
program controllingTurbulence
  implicit none
	
  real*8 q
  real*8 k0 
  real*8 nu 
  real*8 dt
  integer*8 m

  parameter (q = 2)
  parameter (m = 22)
  parameter (dt = 1.d-4)
  parameter (nu = 1.d-7)
  parameter (k0 = 2.d0**(-4))

  complex*16, dimension(1:m):: FORCE
  real*8,     dimension(1:m):: K, K2, A, B, C, LT, NLT
  
  call initForce(FORCE, m)
  call precalculate(K, K2, A, B, C, LT, NLT, nu, k0, q, dt, m)
  call shellModels(FORCE, LT, NLT, A, B, C, dt, m)  
end program controllingTurbulence


! ----------------------------
!
! Computes the energy evolution with the given parameters and outputs this in 
! the file 'velocity.txt'.
!
! ----------------------------
subroutine shellModels(FORCE, LT, NLT, A, B, C, dt, m)
  implicit none 

  integer*8 i
  integer*8 j
  integer*8 n
  integer*8 m
  integer*8 iterBTS
  integer*8 iterOutput
  integer*8 numberLastBTS
  integer*8 numberTransientBTS
  real*8 dt
  real*8 dtol
  real*8 time
  real*8 outputTime
  real*8 bigTimeStep
  real*8,     dimension( 1:m  ):: A, B, C, LT, NLT
  complex*16, dimension( 1:m  ):: FORCE, G0, G1
  complex*16, dimension(-1:m+2):: U
  
  parameter (dtol = 1.d-10)
  parameter (outputTime = 0.1d0)
  parameter (numberLastBTS      = 4)
  parameter (numberTransientBTS = 20)
  parameter (bigTimeStep = 250.D0)
  
  iterBTS = IDINT(bigTimeStep / dt / 2.d0 + dtol)
  iterOutput = IDINT(outputTime / dt / 2.d0 + dtol)
  print *, "Number of iterations in a Big Time Step = ",  iterBTS
  
  
  call initU(U, m)
  call G(G0, U, FORCE, A, B, C, m)
  call G(G1, U, FORCE, A, B, C, m)

  do i = 1, numberTransientBTS
    do j = 1, iterBTS
      call slavedAdamsBashforth2(U, G0, G1, FORCE, LT, NLT, A, B, C, m)
    end do
    
    time = i * bigTimeStep
    print *, i, time
  end do

  open (10, file = 'very_long_velocity.txt')
  do i = 1, numberLastBTS * iterBTS / iterOutput
    do j = 1, iterOutput
      call slavedAdamsBashforth2(U, G0, G1, FORCE, LT, NLT, A, B, C, m)
    end do
    
    write (10, *) time + i * outputTime

    do n = 1, m
      write (10, *) REAL(U(n))
      write (10, *) AIMAG(U(n))
    end do
  end do  
  close (10)
  
  return
end subroutine shellModels


! ----------------------------
!
! Initiates the vector U with proper values.
!
! ----------------------------
subroutine initU(U, m)
  implicit none

  integer*8 m
  complex*16, dimension(-1:m+2):: U

  U    = (0.d0,    0.d0)
  U(5) = (0.03d0, -0.065d0)
  U(6) = (0.06d0, -0.09d0)
  U(7) = (0.06d0, -0.09d0)
	
  return
end subroutine initU


! ----------------------------
!
! Initiates the vector FORCE with proper values.
!
! ----------------------------
subroutine initForce(FORCE, m)
  implicit none
  integer*8 m
  complex*16, dimension(1:m):: FORCE
    
  FORCE = (0.d0,0.d0)
  FORCE(4) = 5.d0 * (1.d0,1.d0) * 1.d-3
  
  return
end subroutine initForce


! ----------------------------
!
! Precalculates the vectors K, K2, A, B, C, LT, NLT
! using the arguments nu, k0, q, dt and m.
!
! ----------------------------
subroutine precalculate(K, K2, A, B, C, LT, NLT, nu, k0, q, dt, m)
  implicit none

  integer*8 m
  integer*8 n
  real*8 q
  real*8 k0 
  real*8 nu 
  real*8 dt
  real*8, dimension(1:m):: K, K2, A, B, C, LT, NLT

  do n = 1, m
    K(n)  = k0 * q**n
  end do
  
  K2 = K**2
  A  = K
  B  = -0.5d0 * K / q
  C  = -0.5d0 * K / q**2
  LT = DEXP(-nu * dt * K2)
      
  if (nu == 0) then
    NLT = dt
  else
    NLT = (1.d0 - LT) / (nu * K2)
  end if
  
  return
end subroutine precalculate
